function [S, V, dxi, W_new] = rcnl_dvarcov( ...
    dparams,           ...
    s_jt,              ...
    xlinear,           ...
    xnonlin,           ...
    cdindex,           ...
    dfull,             ...
    IVD,               ...
    W,                 ...
    IVD_alt,           ...
    W_alt,             ...
    fe,                ...
    xi,                ...
    cluster,           ...
    expmval_mat,       ...
    rho_transform,     ...
    derdefined,        ...
    tolerance          ...
)

dxi        = rcnl_dxi( ...
    dparams,           ...
    s_jt,              ...
    xlinear,           ...
    xnonlin,           ...
    cdindex,           ...
    dfull,             ...
    IVD,               ...
    W,                 ...
    IVD_alt,           ...
    W_alt,             ...
    fe,                ...
    expmval_mat,       ...
    rho_transform,     ...
    tolerance          ...
);
dxi(:, setdiff(1:max(derdefined), derdefined)) = 0;

% check parameter is not on bound
rho = dparams(end);
if rho_transform
    rho = rho_transform * (1/(1+exp(-rho)));
end

nfe    = size(fe, 2);
xtilda = xlinear(:, 1:end-nfe) - fe * ((fe' * fe) \ (fe' * xlinear(:, 1:end-nfe)));

% regular estimator
ME_a   = IVD_alt'*xtilda/length(cluster);
ME_b   = IVD_alt'*dxi/length(cluster);
MI_a   = IVD'*xtilda/length(cluster);
MI_b   = IVD'*dxi/length(cluster);

ME_b = ME_b(:,derdefined);
MI_b = MI_b(:,derdefined);

M0 = [ME_a, ME_b; ...
      MI_a, MI_b];
S0 = [ME_a'*W_alt, zeros(size(ME_a,2),size(W,1)); ...
      zeros(size(MI_b,2),size(W_alt,1)), (MI_b - MI_a*inv(ME_a'*W_alt*ME_a)*ME_a'*W_alt*ME_b)'*W];
SE = [-MI_a*inv(ME_a'*W_alt*ME_a)*ME_a'*W_alt eye(size(W))];

% estimate moment variance matrix
ncol        = size(IVD_alt, 2) + size(IVD, 2);
S           = zeros(ncol, ncol);

size_V      = size(xtilda,2) + size(dxi,2);
V           = zeros(size_V);

ix      = [size(xtilda,2) size(xtilda,2) + derdefined];

for m = unique(cluster)'
    xi_m       = xi(cluster == m);
    IV_m       = [IVD_alt(cluster == m,:), IVD(cluster == m, :)];
    S          = S + IV_m' * (xi_m * xi_m') * IV_m;
end

k      = 1+size(derdefined,2)+size(fe,2);
K      = ((length(cluster)-1)/(length(cluster)-k))*(length(unique(cluster))/(length(unique(cluster))-1));
S      = K*S/length(cluster);


if rho>0.9499
    disp(' ')
    disp('WARNING: rho is close to bound. Variance matrix undefined.')
    disp(' ')
end

V(ix,ix)                = inv(S0*M0)*(S0*S*S0')*inv(M0'*S0')/length(cluster);
W_new                   = pinv(SE*S*SE');

end

